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We investigate electrical transport through a single-electron transistor coupled to a nanomechan- 
ical oscillator. Using a combination of a master-equation approach and a numerical Monte Carlo 
method, we calculate the average current and the current noise in the strong-coupling regime, study- 
ing deviations from previously derived analytic results valid in the limit of weak-coupling. After 
generalizing the weak-coupling theory to enable the calculation of higher cumulants of the current, 
we use our numerical approach to study how the third cumulant is affected in the strong-coupling 
regime. In this case, we find an interesting crossover between a weak-coupling transport regime 
where the third cumulant heavily depends on the frequency of the oscillator to one where it be- 
comes practically independent of this parameter. Finally, we study the spectrum of the transport 
noise and show that the two peaks found in the weak-coupling limit merge on increasing the cou- 
pling strength. Our calculation of the frequency-dependence of the noise also allows to describe how 
transport-induced damping of the mechanical oscillations is affected in the strong-coupling regime. 



I. INTRODUCTION 

Nano-electromechanical devices, i.e., nanostructures in 
which electric transport through a device is influenced by 
its mechanical degrees of freedom and vice versa, have 
attracted a lot of interest recentlyi^ On the one hand, 
these devices are promising for applications like sensors 
or ultra-sensitive mass detectors. On the other hand, 
they have opened up new directions in fundamental re- 
search, with projects to cool nanomechanical systems to 
the quantum limitA^ 

The nanomechanical properties of single-electron tran- 
sistors (SETs) are of particular interest in this context. 
The central island of a SET may be allowed to mechani- 
cally move between the two leads, such that electrons can 
tunnel on the island if the island has approached one lead 
and leave it again once it has mechanically moved to the 
other lead. These "shuttles" have been investigated in 
great detaiL 5 ! 6 ! 7 ! 8 ! 9 ! 10 ! 11 ! 12 ! 13 Another possibility to cou- 
ple the electrical and mechanical properties of the device 
is to design the SET such that its capacitive coupling 
to the gate depends on the displacement of a mechan- 
ical oscillator. Thus, mechanical degrees of freedom of 
the system may strongly influence the current-voltage 
characteristics, current noise, and higher cumulants of 
the currenti 14 i 15 i 16 i 17 i 18 i 19 i 20 It has been shown that the 
Coulomb blockade peaks are split for harmonic oscilla- 
tions and are broadened by thermal oscillations. Knowl- 
edge of the SET transport properties therefore allows one 
to determine the characteristics of the oscillator, such as 
its amplitude and frequency. In such systems, electron 
tunneling through the island also has an effect on the 
motion of the oscillator. This back-action leads to fluc- 
tuations in the oscillator position and to damping^! 

Practical implementations of oscillator-coupled SET 
transistors can be realized by combining nanostructured 
silicon resonators with metallic SETs. 22,23 Another possi- 



bility is to build SETs from suspended carbon nanotubes 
that act as quantum dotsi24 Quite recently, mechanical 
oscillations of the nanotube in such a device have been 
directly observed ~ 

In the following, we will consider a SET transistor cou- 
pled to a classical harmonic oscillator. This system has 
already been studied extensively^**^ However, previous 
studies investigated the regime where the coupling be- 
tween the oscillator and the SET is weak and the ques- 
tion what happens when the coupling is increased is still 
of great theoretical interest^ even if this regime is not 
readily accessible in the current generation of experi- 
ments. In this article, we will use a combination of a 
master-equation approach and a numerical Monte Carlo 
procedure to calculate the electrical current and its sec- 
ond and third cumulants and study how they are mod- 
ified by coupling to the oscillator, in the regime where 
the coupling is strong. We will also study the frequency 
dependence of the transport noise. 

The paper is organized as follows: in Section [HI we 
discuss the system and the model whose strong-coupling 
limit will be studied in the subsequent sections. The 
model and the master-equation approach that we use 
follow closely Ref. Q. This section also introduces the 
important dimcnsionless coupling parameter n that is 
the ratio of the typical mechanical energy scale and the 
source-drain voltage. In the next section, Sec. II I II we 
calculate the probability distributions of the position of 
the oscillator if the SET is in state N or N + 1 using 
a numerical Monte Carlo procedure. The Gaussian form 
predicted by the weak-coupling approach is modified dra- 
matically in the strong-coupling regime. In Section IIVI 
we calculate the average current through the device at the 
charge-degeneracy point of the SET at which the current 
is maximal and discuss the deviations from the weak- 
coupling results. Finally, in Sections IVl and IVTI we extend 
our studies to the current noise and the third cumulant 
of the current. 
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FIG. 1: Circuit diagram of the system studied. The gate 
capacitance of the SET depends on the displacement of a 
mechanical oscillator, leading to a coupling of the electrical 
transport through the device and the mechanical motion of 
the oscillator. 



II. COUPLED SET-OSCILLATOR SYSTEM 
DESCRIPTION 

To describe the coupled SET-nanomechanical oscilla- 
tor system, we use the formalism introduced in Ref. 

The system that we consider is shown in Fig. ^ m 
a schematic way and consists of two symmetric tunnel 
junctions, each with resistance R and capacitance C, con- 
nected in series. Transport through the SET is described 
using the orthodox model, where only two charge states 
are considered and where the current arises only from 
sequential tunneling i2L2& In this case, transport is gov- 
erned by four tunneling rates where i = R, L is the 
lead index and a = +, — indicates the direction of the 
tunneling. In this work, we adopt the convention that 
the forward (+) direction, given by the polarity of the 
bias voltage, is from the right to the left lead. The tun- 
neling rates can be calculated using Fermi's golden rule 
and are a function of the difference in free energy AE of 
the system before and after a tunneling event 



Tf = ^AEf f(AEf) 



(1) 



where f(x) — (1 — e x / fc -B T ) 1 ) with T the electronic 
temperature. The energy differences AEf are given by 



AE+ = -AE7 = eV r , 



' ds (l + (2N-2N g + l)-^_) 



AE^ = -AE R = e 



where Vds is the applied drain-source voltage, E c = 
e 2 /(2C + C g ) is the charging energy of the island and 
N g = CgVg/e is the optimal number of charges on the 
island. Knowing the different rates, the average current 
/ flowing through the SET can be calculated using 



where Pn(n+i) is the probability to find the island in 
charge state N(N + 1) in the stationary limit. 

Our model of the SET remains valid as long as its 
charging energy E c is large compared to the electronic 
thermal energy fcgT e and the source-drain bias eVds- We 
will neglect second-order tunneling processes (cotunncl- 
ing). 

In this work, the nanomechanical oscillator is mod- 
eled as a single, classical, harmonic oscillator of frequency 
luq. Introducing a time scale r t = Re /Vds which has the 
physical meaning of an average time between tunneling 
events, we can use the dimensionless parameter 

e = w T t = w — ( 4 ) 

Vds 

to compare the typical electrical and mechanical 
timescales. 

A particular state of the oscillator is then represented 
by a position x and velocity u. We choose x = to be the 
equilibrium point of the oscillator when iV charges are on 
the SET. When the charge state of the island is changed, 
for example, from N to N + 1, the change in the elec- 
trostatic forces between the oscillator (kept at constant 
potential V g ) and the SET effectively shifts the equilib- 
rium position of the resonator. The distance between the 
equilibrium positions when N and N + 1 charges are on 
the island defines a natural length-scale xq of the prob- 
lem, xo = —2E c N g I (mwgd) . Here, d is the distance sep- 
arating the oscillator's equilibrium position and the SET 
island, such that the gate capacitance depends on x like 
C g (x) ~ (d + x)^ 1 ~ 1 — x/d. From now on, we will also 
use dimensionless rates, i.e., all the rates will be given in 
units of Tj - . 

Coupling a SET and a nanomechanical oscillator sys- 
tem is readily done by using the oscillator itself as the 
SET's gate. In this configuration, the capacitive cou- 
pling between the oscillator and the SET depends on the 
distance between them and by extension on the oscil- 
lator's position, effectively allowing one to monitor the 
dynamics of the oscillator via the SET. As long as the 
amplitude of the oscillations around its equilibrium po- 
sition is small compared to the distance d separating the 
oscillator and the SET island, the gate capacitance C g (x) 
can be treated as linear in i. As a consequence, we ob- 
tain position-dependent dimensionless tunneling rates of 
the form 
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A L = \ + (2N-2N g + l)-^- 
2 eV ds 

A R = \-{2N- 2N g + 1)^- + kN 
2 eVds 



(5) 



(6) 
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are the position-independent part of the full dimension- 
less rate T?(x) that fulfill A L + A R = 1, and 



K = mu Q x /(eV ds ) 



(7) 



is a dimensionless coupling parameter that will play an 
important role in the following. Note that Al, Ar can 
become negative in the strong-coupling limit. The av- 
erage dimensionless current in the presence of position- 
dependent rates can be calculated as an average of the 
different rates weighted by the probability to find the 
oscillator at a position x: 



1= I dx(P N+1 (x)r+{x)-P N (x)Tl(x)) 

) 

dx (P N {x)T+(x) - P N+1 (x)T R (x)) , (8) 



with fjv(jV+i) ( x ) the probability to find the oscillator at 
position x while the island charge state is N(N + 1). 

In the zero-temperature limit, the Fermi functions in 
Eqs. JSJ are in fact Heaviside step functions that deter- 
mine the possible transport direction as a function of the 
position of the oscillator. Indeed, at zero temperature, 
x L = A^xq/k and x R — — ArXq/k define points where 
the current direction at lead L and R changes sign. For 
x R < x current in the right junction can only be di- 
rected towards the island while in the opposite case only 
charge transfer from the island to the right lead is pos- 
sible. Equivalently, transfer through the left junction is 
allowed from the island to the lead if x < x L and from 
the lead to the island otherwise. It is interesting to note 
that transport can be blocked altogether via this mech- 
anism. For example, if N + 1 electrons are on the island 
and the oscillator is in position x > x L , transport of the 
extra charge from the island to any lead is effectively for- 
bidden, our choice of bias direction imposing x R < x L . 

The canonical way of dealing with an SET in the se- 
quential tunneling regime is to introduce a master equa- 
tion for the different charge states of the island. If the 
oscillator is coupled to a nanomechanical oscillator, such 
a simple master equation cannot be written, since the 
tunneling rates depend on the stochastic evolution of the 
oscillator. Following Ref. ^|we can introduce the prob- 
ability distributions Pjsi{x,u;i) and Pn+i {%■> u\ t) to find 
at a time t, the oscillator at position x, u in phase space 
and the SET in charge state N and N + 1 respectively 
and derive a master equation for these new objects: 



d d d 

—P N {x,u;t) = v x—P N (x,u;t) - u—P N (x,u;t) 

+ [T+(x)+T-(x)]P N+ i(x,u;t) 

- [r+(x)+Tl(x)]P N (x,u;t) , (9a) 



d d 
—P N+l (x,u;t) = ujI{x - x ) — P N+ i(x,u;t) 

d 

- u—P N+ i(x,u;t) 

ox 

- [r+^ + r-^)] p N+1 (x,u-,t) 

+ [T+(x)+Tl(x)]P N (x,u;t) . (9b) 

As pointed out in Ref. 0, when the coupling between 
the oscillator and the SET is weak (k <C 1) and when 
the gate voltage V g is such that the system is tuned far 
from the Coulomb-blockade region, one can make the ap- 
proximation that x L — > oo and x R — ► — oo and then write 
the tunneling rates as 



T+(x) = A l ~k- 



r+(x) = a r + k . 



x 

X 



rz(aO = o, 

T-(x) = . 



(10) 



This weak-coupling form of Eq. (JSJ) effectively corre- 
sponds to neglecting any back-currents and the possibil- 
ity of position-induced current blockade. However, the 
master equation is then simple enough to allow analyti- 
cal solutions. 

In this work, we will not study the effect of extrinsic 
damping (i.e., a finite quality factor of the oscillator) and 
of finite temperatures, since they were discussed exten- 
sively in Refs. [TIIT5I 



III. DYNAMICS OF THE OSCILLATOR IN THE 
STRONG-COUPLING REGIME 

In the weak-coupling limit k <C 1, it was founcU* that 
the interaction between the SET and the oscillator in- 
troduces an intrinsic damping mechanism. The damp- 
ing, characterized by a decay rate 7 = ne 2 (in units of 
TjT 1 ) leads to a steady-state solution for the probability 
distributions Pn(n+i)(x, u). In particular, the probabil- 
ity distributions Pn(n+i)( x ) = J du Pn(n+i)( x j m )> from 
which one can calculate the average current, have been 
shown to be well approximated by Gaussians centered at 
x = and x = xq for Pjv and Pn+i-, respectively. 

One of the main goals of this work is to study de- 
viations from the weak-coupling behavior. Without the 
simplifications possible for n -C 1 leading to Eq. (|1U|> , the 
stationary probability distributions Pn(n+i)(x, u) can no 
longer be investigated analytically and numerical meth- 
ods must be used. In this work, we used a Monte 
Carlo approach to simulate the stochastic nature of the 
SET-nanomechanical oscillator system in the parameter 
regime where the typical mechanical energy towqXq is 
comparable to the bias energy eVd s . Details of our im- 
plementation of the Monte Carlo method are given in 
Appendix 1X1 

We study the probability distribution of the oscillator 
in the charge-degenerate case, where (Pn) = (Pn+i), 
where (Pn) — J dxP]^{x). At this point the current 
flowing through the SET is maximal. In the presence of 




FIG. 2: Upper panel: Total probability distribution P(x) = 
Pn(x) + Pn+i(x) of the oscillator's position for different val- 
ues of the coupling constant re (defined in Eq. Q). Lower 
panel: Probability distribution Pn(x) to find the oscillator at 
position x if the SET is in charge state N. Pn+i{x) can be 
obtained by the symmetry relation Pjv+i(a;) = Pjv(1/2 — x). 
In both panels, lines are shifted for clarity by 2k, and the 
difference between neighboring curves is Are = 0.05. All cal- 
culations were done at e = 0.3 and at the charge-degeneracy 
point. For the definition of e, see Eq. 



the oscillator, charge degeneracy is reached when = 
l/2+re/2. This relation, exact in the weak-coupling limit, 
has been empirically verified over the whole range of re 
studied. In the weak-coupling limit, this relation can be 
shown using (Pjv+i) = familiar from 2-state SETs. 
In our case, at the degeneracy point symmetry consider- 
ations impose (x) = xq/2 and the position dependence 
of the rate r^(x) must be taken into account, such that 
(Pjv+i) = I = + k(x/x ). This effectively corre- 
sponds to (Pjv+i) = A L — re/2. 

As can be seen from Fig. |21 as re is increased, the sta- 
tionary position probability distribution evolves continu- 
ously from the weak-coupling Gaussian form to a distri- 
bution showing two sharp peaks at x = and x = x a 
in the limit where re = 1. This evolution is the re- 
sult of a sharpening of each of the two subdistributions 



FIG. 3: Two-dimensional phase-space distributions P(x, u) = 
Pn(x,u) + Pjsr+i(x,u), for k = 0.2 (left panel) and n — 0.9 
(right panel). 



Pn(n+i)(%) around their equilibrium position when re is 
increased, allowing one to resolve the two subdistribu- 
tions individually. This should not only be seen as nat- 
ural consequence of the fact that the typical distance xq 
scales like y^re. In fact, the main cause of the appearance 
of the two sharp peaks is that small amplitude oscilla- 
tions about each of the two equilibrium points become 
very stable when re is increased. 

We also note that the qualitative shape of each sub- 
distribution evolve when re is increased. While at low 
coupling the subdistribution Pjy(x) (resp. Pj^+i(x)) is 
symmetric about x = (resp. x — xq), this is not the 
case for re > 0.4. This asymmetry arises only at higher 
coupling since for low re, the probability to find the oscil- 
lator at x < x R or x > x L is negligible. When re > 0.4, 
the probability of the oscillator being in a region trans- 
port is forbidden becomes important. Symmetry break- 
ing arises since these regions are located only on one side 
of each equilibrium point. 

Finally, we note that the important changes in 
Pn(n+i)( x ) that accompany a variation in re are 
also seen in the stationary velocity subdistributions 
Pv(jv+i)( u ) = J <^' u Pv(iv+i)(a ; , u )j that approximatively 
follow P/v(u/euo) — Pn(x/xq) and P/v+i(w/ eu o) — 
Pn+i((x — xq)/xq), where uq — xo/r t is the typical ve- 
locity scale in the problem. This can be seen in the two- 
dimensional phase-space distributions shown in Fig. |3J 
where both P(u) and P(x) are shown to become more 
peaked when re is increased. 

For re > 1, our numerical investigations show that the 
current is strongly suppressed, rendering the intrinsic 
damping mechanism discussed at the beginning of this 
section ineffective. In this case, the system cannot be 
characterized by a steady-state probability distribution, 
and our model is not appropriate. Therefore, we only 
studied the parameter range re < 1. 

A similar reasoning applies to Coulomb-blockade re- 
gion, where the damping of the oscillator's motion is 
severely suppressed. However, numerically finding a 
steady-state solution close to the degeneracy point is pos- 
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FIG. 4: Current (in units of e/rt) as a function of k at the de- 
generacy point (Pjv) = (Pjv+i), for e = 0.3. The dots are the 
results of the Monte Carlo calculation and the solid line is the 
analytic form found within the weak-coupling approximation. 



sible. 



IV. AVERAGE CURRENT 

The average current flowing through the SET is closely 
tied to the oscillator's position distributions via the 
position-dependent tunneling rates. Consequently, one 
expects that the deviations from the weak-coupling be- 
havior observed in P(x) = Pn(x) + Pn+i{x) would affect 
the current characteristics when n is increased. 

Just like in the previous section, we focus on the de- 
generacy point where the average charge state of the is- 
land is N + 1/2. At this point, the weak-coupling the- 
ory predicts^ that the current decreases linearly with in- 
creasing k: I — i^(e/r t ). This decrease in the current 
can be explained in a qualitative way by the reduction 
of the overlap of the distributions Pn(x) and Pn+i(x) 
as the coupling is increased, each distribution becoming 
more localized around its equilibrium point, see Fig. [21 

Figure 21 shows the average current as a function of 
k. Like in the weak-coupling limit, the localization of 
the different probability distributions around its equilib- 
rium points leads to an overall decrease in the current 
when the coupling grows stronger. For k > 0.3, however, 
we see that the numerical results deviate from the weak- 
coupling behavior: for stronger coupling the current is 
higher than the weak-coupling result. This can be ex- 
plained using the rates given by Eq. (|1 . When either 

dx Pn+i or J^^dx Pn is not negligible, these rates 
allow unphysical backward currents that are not present 
in the full master equation. For example, a point located 
at x > x L in the steady-state probability distribution 
Pn+i{x,u) would contribute negatively to the average 
current when using the rates calculated within the weak- 
coupling approximation while it would not contribute to 
the current when taking into account the full expression 



for the rates given in Eq. ©. 

Over the range of frequencies that we studied numeri- 
cally (0.1 ^ e ^ 0.4), the current was found to be practi- 
cally independent of e for all but the strongest couplings 
(k > 0.8). For instance, at k = 0.9, the difference be- 
tween the calculated currents at e = 0.1 and e = 0.4 is 
on the order of a few percent. 



V. NOISE AND HIGHER CUMULANTS IN THE 
STATIC REGIME 

Originally, the interest in SETs was motivated by 
the suppression of the current in the Coulomb-blockade 
regime and the high sensitivity of the current to small 
variations of the gate voltage. However, it is clear that a 
complete description of the transport processes in these 
devices requires not only knowledge of the current, but 
also of current-current correlations like e.g., the cur- 
rent noisei^ii^i Recently, higher-order correlations have 
also been studied both theoretically and experimentally 
in nanoscale devices, in the framework of full counting 
statistics (FCS) (see R ef. |32l for a collection of articles on 
this topic, and Ref. I3ll3l for a description of FCS in the 
context of transport through SETs). The FCS approach 
consists in studying the probability distribution V n (to) 
that n electrons are transferred through one lead of the 
SET within a time period to, in the limit where to is by 
far the longest time scale in the problem. The full infor- 
mation about the transport properties is contained in the 
cumulants of this distribution function, the first three of 
which are given by the average \x\ = (n), the variance 
= (n 2 ) — ( n ) 2 , and the skewness ^3 = ((n — (n)) 3 ) 
that measures the asymmetry of the distribution. For 
example, the current / = e(n) /to is proportional to the 
mean of this distribution, while the zero-frequency shot 
noise power S(0) = 2e 2 /Z2/to is determined by its second 
cumulant. 

In this section, we study in detail the second and third 
cumulants of the probability distribution function V n (to) 
in the case of a coupled SET-nanomechanical system. 



A. Weak-coupling case 

It is instructive to start by considering the weak- 
coupling case k < 1, since in this regime we can cal- 
culate the noise and higher cumulants without resort- 
ing to Monte Carlo simulations by solving directly for 
(n l (io)) in the long-time limit (toj ^> 1). In this section, 
we generalize the work that was done in Ref. where a 
method to calculate the current-noise using the moments 
of the steady-state probability distribution Pn(n+i) ( x > u ) 
of the oscillator in phase space was described. In this ap- 
proach, the current-noise is calculated from the solution 
of the equation of motion of (n 2 (£)), the average square 
of the number of charges that went through a lead in a 
time t. Here, we extend this method for the calculation 
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of higher cumulants by deriving the equation of motion 
for the general quantity (n m (t)} from which the m— th 
cumulant can be extracted. 

To proceed, we write a master equation for the prob- 
ability -PjvYjv+i) ( x i u i to find, at time t, the oscillator 
at position x with velocity u, the island being in charge 
state N(N+1), and that n charges having passed through 
a lead of the SET in the interval [0;t]. We will again 
make the assumptions leading to Eq. I jlOj l. Considering 
for definitiveness the left lead, at zero temperature and 
neglecting any extrinsic damping, one findsi 5 - 



d d d 

— P%(x, u; t) = ufa—P^x, u; t) - U —P£(x, u; t) 

+ T+(x)Pfc\(x, u; t) - T+(x)P^(x, u- 1) , 

(11a) 

d d 
T^P r N +1 {x,u;t) = uj 2 (x - x )—P^ +1 (x,u;t) 

d 

- Ug-P N+1 (x,u;t) 

- T+ (x)P% +1 (x, u: t) + T+ (x)P% (x, u; t) . 

(lib) 

where the rates are taken from Eq. (|10|l . Defining the 
coupled moments (x 3 u k n m ) and (x 3 u k n m ) ^r+i as 

(n m x 3 u k ) = 

^n m du dxx j u k [P%(x,u;t) + P^ +1 (x,u;t)] , 

n ■* 

(12a) 

{n m x j u k ) N+ -L = y^n m I du I dxx 3 u k P^ +1 (x,u;T) , 

n J ^ 

(12b) 

one can calculate the equation of motion for these quan- 
tities using Eq. I|ll(l . With x in units of xq and u in units 
of uq, one finds 

^(x 3 u k n m ) N+1 = 
dr 

- fee 2 [{x J+1 u k - 1 n m ) N+1 - (.T^ fe - 1 n m ) w+1 ] 
+ j(x 3 - 1 u k+l n m ) N+1 - (x 3 u k n m ) N+1 

+ A R (x 3 u k n m ) + n(x 3+l u k n m ) , (13a) 



(x 3 u k n m ) - 



dr 



- ke 2 [(x 3+l u k ~ l n m ) - (x j u k - 1 n m ) N+1 ] 
+ j(x 3 ^ 1 u k+1 n m ) 

m—l , v 

i m ) [Ai{i 3 'u l nV +1 -K(^ 1 u*n% +1 ] 



Here, averages that are independent of n (averages of the 
form (x 3 u k n )) are time-independent and can be evalu- 
ated in the stationary limit, i.e., Eqs. H13all3b(l can be 
used to generate a closed linear system of equations. 35 
The terms (x 3 u k ) of order j + k — c are connected to the 
terms (x 3 u ) m+i of order j + k = c— 1. This means that 
to solve for a moment (x 3 u k ), we must use the c+ 1 equa- 
tions of the type of Eq. (jl3b(l where j + k = c and the c 
equations of the type of Eq. I|13a|) where j + k = c — 1. 
This method can be used to calculate any moment of 
the form (x 3 u k ) and (x 3 u k )M+i- Knowledge of (x 3 u k ) 
enables one to calculate the long-time behavior of the 
coupled moments of the charge and oscillator's position 
in phase space (x 3 u k n l ), and thus the i-th moment (n l ) 
of TV 

The ratio of the zero-frequency shot noise power and 
the average current (times 2e) , or equivalently, the ratio 
of the second and first cumulants of V n is called Fano fac- 
tor and is readily calculated using this approach. Since 
it shows a complex dependence on the coefficients Al, 
A#, and on the parameters k and e, it is convenient to 
expand the result in powers of n. Introducing a param- 
eter a defined via a = Al — (1 + n)/2 (or, equivalently, 
a = — A R + (1 — k)/2) that measures the difference be- 
tween Al and its value at the degeneracy point, one can 
write down the Fano factor in a way that underlines its 
symmetry with respect to this point: 
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For e <C 1, the Fano factor is dominated by the term 
~ k/e 2 , like in the case where one considers a system 
of two SETs coupled by an oscillator i2i Finally, we note 
that current conservation implies that the Fano factor is 
identical in both leadsi^I 

Equation 1131 is one of the main result of our article, 
as it allows for the calculation of higher cumulants of 
the current by integrating the equation of motion for the 
moments of the form (x 3 u k n m (t)) with m > 0. For exam- 
ple, we calculated the normalized third cumulant [13/(11) 
of V n (to)- The results are presented in Fig.|SJ We stress 
that these results have been obtained by integrating the 
equation motion for (n 3 (t)} valid in the weak-coupling 
regime and not via a Monte Carlo simulation. Starting 
from the value 1/4 at k — O^ 3 . the normalized third cu- 
mulant decreases rapidly when k is increased. On further 
increase of k, it reaches a minimum at an e-dependent 
value of k. The inset of Fig.JSJshows that the leading con- 
tribution to the normalized third cumulant in the weak- 
coupling regime is of the form e~ 4 . As a consequence, we 
note that the asymmetry of the probability distribution 
V n that is determined by /13, can effectively be tuned 
by changing the frequency of the oscillator or r t . The 
scaled quantity e 4 ^3/(71) shows contributions of higher- 
order corrections in e to the normalized third cumulant. 
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FIG. 5: Normalized third cumulant as a function of k for 
different values of e, as calculated within the weak-coupling 
approximation and scaled by e 4 . Dotted line: e = 0.4, dashed 
line: e = 0.3, solid line: e = 0.2, dash-dotted line: e = 0.1. 
The inset shows the e-dependence of the normalized third 
cumulant at k = 0.1 (symbols), as calculated within the weak- 
coupling approximation. The solid line is a fit to a power 
law ~ e~ 4 . These results were obtained by integrating the 
equation of motion for (n 3 (i)} following from Eq. 1131 . 



B. Higher coupling 

It is unfortunately not straightforward to generalize 
the previously described method for calculating the cu- 
mulants of V n outside the weak-coupling regime. The 
presence of x-dependent Fermi functions in the tunnel- 
ing rates as well as the possibility of charge flow against 
the direction set by the bias voltage due to the position 
of the oscillator gives rise to a system of equations that is 
not closed and cannot be solved analytically. Even if we 
neglect transport against the dominant direction of the 
current Tj(x) ~ 0, but keep the position dependence of 
the Fermi distributions in r+(x), it is still not possible 
to derive a system of equations coupling only objects of 
the form (x^u ri m ). Therefore, we will use a numerical 
approach to evaluate the cumulants of V n . 

Indeed, the Monte Carlo method described in Ap- 
pendix [S] can be used to measure the FCS of elec- 
tron transport in the same way as it can be done in 
experiments. 38 A very long Monte Carlo run is divided 
into intervals of duration to 7 , here, 7 = ne 2 is the 
damping constant defined at the beginning of Sec. Ill II 
7 _1 is the longest intrinsic time scale of the problem. By 
counting the number of charges going through one lead 
during each interval, one can reconstruct the probability 
distribution V n (to), and from it calculate the cumulants. 

We study current correlations at the charge degener- 
acy point, where the average charge state of the island 
is N + 1/2. The top panel of Fig. compares the weak- 
coupling Fano factor to the numerical Monte Carlo re- 
sults for different values of the coupling parameter k. 
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FIG. 6: Upper panel: Fano factor as a function of n at the 
degeneracy point (Pjv) = (Pjv+i). The dots are the results 
of the numerical calculation and the solid line is the analytic 
form found within the weak-coupling approximation. Lower 
panel: Normalized third cumulant ^13 of the probability dis- 
tribution V n . For both panels, e = 0.3. 



Naturally, for k < 0.2, the agreement between the numer- 
ical results and those obtained analytically is very good. 
Beyond this point, the numerically calculated Fano fac- 
tor shows an interesting non-monotonic behavior, with 
a maximum at k ~ 0.35 and a minimum at k ~ 0.85. 
The lower panel of Fig. also shows the evolution of the 
normalized third cumulant of V n , giving the asymme- 
try of this probability distribution about its mean (n). 
Starting from the k = value of 1/4 derived for a sim- 
ple SET device, our results show that this quantity is, 
in the weak-coupling limit, very sensitive to variations of 
k. Indeed, the normalized third cumulant changes sign 
twice in the region k < 0.35, reaching a maximum value 
approximatively in the middle of this region. This con- 
trasts with the strong-coupling behavior: [13/(71) stays 
practically constant for 0.5 < n < 0.9. 

We will now address the question how the previous re- 
sults are modified when changing the frequency of the 
oscillator. Figure [7| shows the dependence of the Fano 
factor and the normalized third cumulant as a function 
of k for different values of e. First, we note that the 
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FIG. 7: Fano factor (upper panel) and normalized third cu- 
mulant (lower panel and inset) as a function of k for different 
values of e: e = 0.1 (squares in the upper panel and in the 
inset), e = 0.2 (diamonds), e = 0.3 (circles), and e = 0.4 
(triangles). Note the logarithmic y-axis in the upper panel. 



actual value of the Fano factor is increased dramatically 
for lower oscillator frequencies, as expected from the term 
~ n/e 2 that dominates in the low- frequency regime. In 
the weak-coupling region (re < 0.3), the magnitude of the 
normalized third cumulant is also heavily affected by a 
change in frequency, in agreement with the weak-coupling 
leading order dependence ~ e~ 4 . Despite these major 
changes in magnitude of both the Fano factor and the 
normalized third cumulant, the overall qualitative effect 
of an increase in coupling does not seem to depend heav- 
ily on e, in the frequency range we investigated. In par- 
ticular, the position of the maximum in the Fano factor 
remains constant. Also, the normalized third cumulant 
always shows a change of sign, albeit at an e-dependent 
value of re, and goes to a positive for re — > 1. Remarkably, 
the value of the normalized third cumulant is much less 
sensitive to e in the strong-coupling regime. This might 
be the signature of a transport regime in the re — > 1 region 
that is radically different of the one found for re ~ 0.2. 



VI. FREQUENCY-DEPENDENT NOISE 

In systems that exhibit internal dynamics like the 
one we study, it is especially interesting to look at 
the frequency-dependence of the current-current corre- 
lations. In Ref. 15, the frequency-dependent noise S(lo) 
of a SET weakly coupled to a nanomechanical oscilla- 
tor was thoroughly studied. It was found that, the noise 
spectrum shows only two peaks at finite frequency at 
oj' and 2u>' , where loq = o>o\/l ~~ K is the effective fre- 
quency of the damped harmonic oscillator. In this sec- 
tion, we extend the work of Ref. by studying the 
frequency-dependent noise power S(u>) in the strong- 
coupling regime (0.2 < re < 1). 

In order to calculate the frequency-dependent noise us- 
ing our Monte Carlo method, we follow the approach 
developed by MacDonald^i*^ that was used recently to 
study the noise properties of mesoscopic systems, includ- 
ing coupled SET-nanomechanical systems in the weak- 
coupling regime. In general, the current-noise power at 
junction a is defined as the Fourier transform of the cur- 
rent autocorrelation function K% t i at junction i, 



cLt cos(wr)ifi.i(r) 



where 



K iti (T) = (I i (t + T)I i (t))-(I i )< 



(15) 



(16) 



To proceed with the MacDonald approach, Jj, and there- 
fore li — (ij), must be assumed to be statistically fluc- 
tuating variables, such that the autocorrelation function 
Ka is independent of t. In this case, the MacDonald 
formula relates the fluctuation Sn about the average of 
the number of charges n that went through a junction in 
time r, 

5m(T) = n t {r) - (h)T = f + dt' {1^) - {h)) , (17) 



to the current-noise power via 







S iti (u) = 2loJ dr sin^r) ^-((6n t ( T )) 2 ) , (18) 

where ((Sn^r)) 2 } = (n 2 (r)) - (l t ) 2 T 2 . Since (n 2 (r)) and 
(Jj) are easily accessible through the Monte Carlo simula- 
tion, S{ui) can be calculated by taking a numerical time- 
derivative of (((5nj(r)) 2 ) and then evaluating the Fourier 
sine-transform. Note that we consider only the particle 
current fluctuations here. The electrical current noise at 
finite frequencies includes a contribution from displace- 
ment currents, which depend on the capacitive couplings 
between the island and the leads^ Since we assume that 
our frequencies of interest are much smaller than the rel- 
evant RC-frequencies, we can neglect the displac ement 
currents here, see e.g. the discussion in Ref. l40l4lL 

The results of the Monte Carlo simulation are shown 
in Fig. |SJ Like in the weak-coupling case, S(u>) shows 
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FIG. 8: Frequency-dependent noise power beyond the weak- 
coupling approximation. For each curve, the SET is tuned to 
the charge-degeneracy point and e = 0.3. 



two main finite- frequency features. Surprisingly, even for 
strong coupling, we do not find any features for frequen- 
cies higher than 2u>o. We find a low- frequency peak that 
defines the frequency co' (which will be different from the 
weak-coupling prediction tooy/1 — k in the general case). 
The second feature evolves from the peak located at 2ui' 
predicted by the weak-coupling theory. While both peaks 
are considerably broadened by an increase in the coupling 
strength, their respective shapes evolve in a qualitatively 
different way. Whereas the first peak is shifted in abso- 
lute frequency, the second peak broadens in a very asym- 
metric way, with much of the weight shifting to lower fre- 
quencies. The slope of its left shoulder decreases with in- 
creasing k until it forms a plateau at around k ~ 0.7. On 
increasing k even further, the two peaks merge, leading to 
super-Poissonian frequency-dependent noise throughout 
the frequency range ui < 1.5wo- 

Figure shows the position of the maxima of the first 
peak in the frequency-dependent noise power as a func- 
tion of k. By comparing the position of the first peak 
extracted from the curves shown in Fig. [S] (data points in 
Fig. to the weak-coupling prediction lu'q = — k 

(solid line in Fig.[§J), we find quantitative agreement only 
for k < 0.2. Beyond this point, the ratio oj' /ojq still de- 
creases, albeit slowly, when k is increased. It reaches a 
saturation value co' ~ 0.7wo for k > 0.7. 

This behavior can be understood by interpreting the 
frequency shift in terms of an effective damping mech- 
anism caused by electron tunneling. Since there is no 
damping without current, the natural modification of the 
weak-coupling damping constant 7 = ne 2 in the strong- 
coupling regime is to renormalize the weak-coupling re- 
sult by the probability P* to find the oscillator in a posi- 
tion where in principle current is allowed, i.e., for x R < x 
and charge state N, or x < x L and charge state N + 1. 
Defining a renormalized damping constant r ) sc = P*ne 2 , 
it is possible to estimate the position of the first peak as 
a function of k using values of P* extracted from curves 
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FIG. 9: Position ui' of the first peak in the frequency- 
dependent noise power as a function of k. The solid line gives 
the weak-coupling prediction Uos/1 — k, the data points are 
the numerical Monte Carlo results, and the dashed line was 
obtained using an effective damping constant, see text. 



presented in Fig. |2 The result is shown as the dashed 
curve in Fig. |5| and agrees with the data points obtained 
by the Monte Carlo method in a quantitative way. 



VII. CONCLUSION 

In this paper, we have studied the strong-coupling limit 
of a SET transistor coupled to a classical harmonic oscil- 
lator. We have used a combination of a master-equation 
approach and a numerical Monte Carlo procedure to 
calculate the position distribution of the oscillator, the 
electrical current, and the zero- frequency noise in both 
the weak-coupling and strong-coupling regime. With in- 
creasing coupling, we found that the position distribu- 
tion of the oscillator evolves from a broad Gaussian to 
a a function sharply peaked around each of the charge- 
state dependent equilibrium positions of the oscillator. 
We found that the average current in the strong-coupling 
regime is higher than the value predicted by the weak- 
coupling theory and that the Fano factor varies in a non- 
monotonous fashion when coupling is increased. We have 
generalized the weak-coupling theory to allow the cal- 
culation of higher cumulants of the current, and have 
presented results for the third cumulant. In the weak- 
coupling regime, the third cumulant was found to de- 
pend strongly on the frequency of the oscillator, whereas 
in the strong-coupling regime it becomes practically in- 
dependent of this parameter. We have also studied the 
frequency-dependent transport noise. Even in the strong- 
coupling regime, there are no peaks for frequencies higher 
than 2wo, and the two peaks found in the weak-coupling 
limit merge on increasing the coupling strength. Fi- 
nally, we introduced a generalized expression connecting 
the damping rate in the strong-coupling regime with the 
other parameters of our model and used it to understand 
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the evolution of the oscillator's damping-renormalized 
frequency as a function of coupling. 
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APPENDIX A: DETAILS OF THE MONTE 
CARLO APPROACH USED 

Monte Carlo methods have been used for a long time 
to calculate numerically the transport properties of meso- 
scopic systems like SETs42, When dealing with a simple 
SET system, the idea of the Monte Carlo approach is 
to solve the master equation for the charge states of the 
SET by discretizing time into small intervals and allow- 
ing charge transfer to and from the dot with a probability 
that is proportional to the tunneling rates and the time 
interval between two attempts. 

If the SET is coupled to a harmonic oscillator, we can 
proceed in a similar way, by considering charge transfer 
attempts at a finite frequency (Art)" 1 , where A -C 1 
is a dimensionless step size. The success probability 



for a charge transfer is calculated from the oscillator's 
position-dependent instantaneous rates T?(x) calculated 
at the time of the attempt. Between different attempts, 
the oscillator's state evolves according to the classical 
equation of motion, whose solution depends on the charge 
state of the SET, the equilibrium position of the oscilla- 
tor being shifted by xq when the charge state is changed 
from TV — > TV + 1 or by — xq in the opposite case. At the 
beginning of each run, the state of the system is deter- 
mined randomly from the stationary probability distri- 
butions Pn{x,u) and P^ +1 {x,u). In practice, this can 
be implemented by using the final state of the (n — l)-th 
Monte Carlo run as the initial state of the n-th run. 

We consider runs of total time ton, such that each run 
consists of M = to/A steps. Both the time scales ton 
and Ar t are chosen in a way such that increasing t or 
decreasing A does not affect the value of the different 
physical quantities we extract from our calculation. In 
practice, this corresponds to choosing A < 0.02 and t n 
an order of magnitude greater than the typical damping 
time I/7. A consequence of this last constraint is that 
the Monte Carlo approach is particularly useful in the 
strong-coupling regime, where the number of steps M 
per run can be kept relatively small, allowing for more 
runs to be made in the same amount of computer time. 

We checked that our code reproduces the analytical re- 
sults of Ref. for the dependence of the noise and the 
third cumulant on the asymmetry coefficients A^ — Ar 
in the K = limit. Also, the probability distributions 
Pn{x,u) and Pn+i{x,u) that we calculate using the 
Monte Carlo approach coincide with the results one finds 
when solving Eq. @ on a gridii^ 
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